QK Box 11.1

Stokes et al. (2014) studied the neurotoxin tetrodotoxin (TTX) in flatworms. The between-plots factor was flatworm species (fixed with two groups: Bipalium adventitium and Bipalium kewense) with individual flatworms (plots) nested within species. The within-plots factor was body segment (fixed with three groups: head, anterior body, posterior body) and each segment represented a “sub-plot”. The response variable was the TTX concentration of tissue adjusted for weight. The main research questions were about the fixed effects of species, body segment and their interaction on TTX concentration, but the analyses also provide information about the variances associated with the random effects of individual within species and the random interaction between individuals within species and body segment.

Bipalium advenitium. Sanjay Acharya, CC BY 4.0, via Wikimedia Commons

Bipalium kewense. Jean-Lou Justine​, Leigh Winsor, Delphine Gey, Pierre Gros, and Jessica Thévenot, CC BY 4.0, via Wikimedia Commons

The data are here

Stokes, A. N., Ducey, P. K., Neuman-Lee, L., Hanifin, C. T., French, S. S., Pfrender, M. E., Brodie, E. D., 3rd & Brodie, E. D., Jr. (2014). Confirmation and distribution of tetrodotoxin for the first time in terrestrial invertebrates: two terrestrial flatworm species (Bipalium adventitium and Bipalium kewense). PLoS One, 9, e100718.

Preliminaries

First, load the required packages (afex, car, lattice, lme4, lmerTest, nlme, VCA, ez, emmeans) and, for convenience, apaTables

Import stokes data file (stokes.csv)

stokes <- read.csv("../data/stokes.csv")
head(stokes,10)

Set contrasts from afex make individual a factor, make sure species a factor too put segments into sensible order

set_sum_contrasts()

stokes$indiv <- factor(stokes$indiv)
stokes$species<-factor(stokes$species)
stokes$segment <- factor(stokes$segment, levels=c("h","b","p"))

Check residuals by leaving out error term

stokes1.aov <- aov(ttxweight~species*segment, stokes)
plot(stokes1.aov)

Wedge-shaped with mean-variance relationship - redo after log transform

stokes$lttxweight <- log10(stokes$ttxweight)
stokes2.aov <- aov(lttxweight~species*segment, stokes)
plot(stokes2.aov)

Fit full model with log(ttxweight)

stokes3.aov <- aov(lttxweight~species*segment+Error(indiv), stokes)
summary(stokes3.aov)

Error: indiv
          Df Sum Sq Mean Sq F value Pr(>F)  
species    1  0.809   0.809    6.17  0.032 *
Residuals 10  1.311   0.131                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
                Df Sum Sq Mean Sq F value  Pr(>F)    
segment          2  17.15    8.57    47.0 2.8e-08 ***
species:segment  2   0.07    0.04     0.2    0.82    
Residuals       20   3.65    0.18                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Use ez for comparison with type 3 SS - same result as design is balanced

ezstokes <- ezANOVA(data=stokes, dv=lttxweight, wid=indiv, within=segment, between=species, type=3)
print(ezstokes)

$ANOVA
           Effect DFn DFd      F        p p<.05    ges
2         species   1  10  6.170 3.23e-02     * 0.1401
3         segment   2  20 46.965 2.78e-08     * 0.7756
4 species:segment   2  20  0.202 8.19e-01       0.0147

$`Mauchly's Test for Sphericity`
           Effect     W     p p<.05
3         segment 0.719 0.227      
4 species:segment 0.719 0.227      

$`Sphericity Corrections`
           Effect   GGe    p[GG] p[GG]<.05   HFe    p[HF] p[HF]<.05
3         segment 0.781 6.94e-07         * 0.899 1.21e-07         *
4 species:segment 0.781 7.65e-01           0.899 7.96e-01          

GG issue, so print clean table with corrections

apa.ezANOVA.table(ezstokes, correction="GG")

ANOVA results
 

Note. df_num indicates degrees of freedom numerator. df_den indicates degrees of freedom denominator. 
Epsilon indicates Greenhouse-Geisser multiplier for degrees of freedom, 
p-values and degrees of freedom in the table incorporate this correction.
ges indicates generalized eta-squared.
 

Get var components using OLS

Note that these estimates treat B(A)*C as the residual for B(A) vc

stokes.vca <- anovaMM(lttxweight~species/(indiv)+segment+species*segment, NegVC=TRUE, stokes)
stokes.vca

ANOVA-Type Estimation of Mixed Model:
--------------------------------------

    [Fixed Effects]

                         int          speciesBadventitium              speciesBkewense                     segmentb 
                     -0.3854                       0.3262                       0.0000                       0.2727 
                    segmenth                     segmentp speciesBadventitium:segmentb     speciesBkewense:segmentb 
                      1.6674                       0.0000                       0.0689                       0.0000 
speciesBadventitium:segmenth     speciesBkewense:segmenth speciesBadventitium:segmentp     speciesBkewense:segmentp 
                     -0.1482                       0.0000                       0.0000                       0.0000 


    [Variance Components]

  Name          DF        SS       MS       VC       %Total     SD       CV[%]     
1 total         29.367096                   0.165407 100        0.406702 98.914014 
2 species:indiv 10        1.310664 0.131066 -0.01717 -10.380622 0        0         
3 error         20        3.651543 0.182577 0.182577 110.380622 0.42729  103.921222

Mean: 0.411 (N = 36) 

Experimental Design: balanced  |  Method: ANOVA

VCAinference(stokes.vca, alpha=0.05, ci.method="satterthwaite")

Inference from Mixed Model Fit
------------------------------

> VCA Result:
-------------

    [Fixed Effects]

                         int          speciesBadventitium              speciesBkewense                     segmentb 
                     -0.3854                       0.3262                       0.0000                       0.2727 
                    segmenth                     segmentp speciesBadventitium:segmentb     speciesBkewense:segmentb 
                      1.6674                       0.0000                       0.0689                       0.0000 
speciesBadventitium:segmenth     speciesBkewense:segmenth speciesBadventitium:segmentp     speciesBkewense:segmentp 
                     -0.1482                       0.0000                       0.0000                       0.0000 


    [Variance Components]

  Name          DF      SS     MS     VC      %Total   SD     CV[%]   
1 total         29.3671               0.1654  100      0.4067 98.914  
2 species:indiv 10      1.3107 0.1311 -0.0172 -10.3806 0      0       
3 error         20      3.6515 0.1826 0.1826  110.3806 0.4273 103.9212

Mean: 0.411 (N = 36) 

Experimental Design: balanced  |  Method: ANOVA


> VC:
-----
              Estimate     DF CI LCL CI UCL One-Sided LCL One-Sided UCL
total           0.1654 29.367  0.105  0.298         0.113         0.270
species:indiv  -0.0172  0.784                                          
error           0.1826 20.000  0.107  0.381         0.116         0.337

> SD:
-----
              Estimate     DF CI LCL CI UCL One-Sided LCL One-Sided UCL
total            0.407 29.367  0.324  0.546         0.336          0.52
species:indiv    0.000  0.784                                          
error            0.427 20.000  0.327  0.617         0.341          0.58

> CV[%]:
--------
              Estimate     DF CI LCL CI UCL One-Sided LCL One-Sided UCL
total             98.9 29.367   78.9    133          81.7           126
species:indiv      0.0  0.784                                          
error            103.9 20.000   79.5    150          82.9           141


95% Confidence Level  |  CIs for negative VCs excluded  
Satterthwaite methodology used for computing CIs 

Fit mixed effects models with lmer

stokes.lmer <- lmer(lttxweight~species+segment+species*segment+(1|indiv), REML=TRUE, stokes)

boundary (singular) fit: see help('isSingular')

Check residuals from lmer plot

plot(stokes.lmer)

summary(stokes.lmer, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: lttxweight ~ species + segment + species * segment + (1 | indiv)
   Data: stokes

REML criterion at convergence: 50.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8114 -0.3996  0.0869  0.7268  1.7277 

Random effects:
 Groups   Name        Variance Std.Dev.
 indiv    (Intercept) 0.000    0.000   
 Residual             0.165    0.407   
Number of obs: 36, groups:  indiv, 12

Fixed effects:
                  Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)         0.4112     0.0678 10.0000    6.07  0.00012 ***
species1            0.1499     0.0678 10.0000    2.21  0.05146 .  
segment1            0.9598     0.0959 20.0000   10.01  3.1e-09 ***
segment2           -0.3264     0.0959 20.0000   -3.40  0.00281 ** 
species1:segment1  -0.0609     0.0959 20.0000   -0.64  0.53257    
species1:segment2   0.0477     0.0959 20.0000    0.50  0.62447    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) specs1 sgmnt1 sgmnt2 spc1:1
species1     0.000                            
segment1     0.000  0.000                     
segment2     0.000  0.000 -0.500              
spcs1:sgmn1  0.000  0.000  0.000  0.000       
spcs1:sgmn2  0.000  0.000  0.000  0.000 -0.500
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

anova(stokes.lmer, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
                Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)    
species           0.81    0.81     1    10    4.88   0.052 .  
segment          17.15    8.57     2    20   51.84 1.2e-08 ***
species:segment   0.07    0.04     2    20    0.22   0.802    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Contrast head vs anterior and anterior vs posterior segments

stokes.emm <- emmeans(stokes.lmer, ~segment)

NOTE: Results may be misleading due to involvement in interactions

stokes.con <- contrast(stokes.emm, "consec")
summary(stokes.con, adjust="none")

 contrast estimate    SE df t.ratio p.value
 b - h      -1.286 0.166 20  -7.750  <.0001
 p - b      -0.307 0.166 20  -1.850  0.0792

Results are averaged over the levels of: species 
Degrees-of-freedom method: kenward-roger 

Get variance components with CIs

stokes.ci <- confint.merMod(stokes.lmer)

Computing profile confidence intervals ...
Warning: Last two rows have identical or NA .zeta values: using minstepWarning: non-monotonic profile for .sig01Warning: bad spline fit for .sig01: falling back to linear interpolationWarning: collapsing to unique 'x' values

stokes.vc <- (stokes.ci)^2
print(stokes.vc)

                    2.5 % 97.5 %
.sig01            0.00000 0.0499
.sigma            0.08976 0.2274
(Intercept)       0.08213 0.2870
species1          0.00064 0.0753
segment1          0.61411 1.2906
segment2          0.25256 0.0226
species1:segment1 0.05620 0.0133
species1:segment2 0.01652 0.0501